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Abstract: The Horizon Run 4 is a cosmological IV-body simulation designed for the study of coupled 
evolution between galaxies and large-scale structures of the Universe, and for the test of galaxy formation 
models. Using 6300'^ gravitating particles in a cubic box of Lbox = 3150 /i“^Mpc, we build a dense 
forest of halo merger trees to trace the halo merger history with a halo mass resolution scale down to 
Mg = 2.7 X We build a set of particle and halo data, which can serve as testbeds for 

comparison of cosmological models and gravitational theories with observations. We find that the FoF 
halo mass function shows a substantial deviation from the universal form with tangible redshift evolution of 
amplitude and shape. At higher redshifts, the amplitude of the mass function is lower, and the functional 
form is shifted toward larger values of ln(l/CT). We also find that the baryonic acoustic oscillation feature 
in the two-point correlation function of mock galaxies becomes broader with a peak position moving to 
smaller scales and the peak amplitude decreasing for increasing directional cosine /i compared to the linear 
predictions. From the halo merger trees built from halo data at 75 redshifts, we measure the half-mass 
epoch of halos and find that less massive halos tend to reach half of their current mass at higher redshifts. 
Simulation outputs including snapshot data, past lightcone space data, and halo merger data are available 
at http://sdss.kias.re.kr/astro/Horizon-Run4/. 
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1. Introduction 

Over the last decade, a series of cosmological N- 
body simulations called the Horizon Run (HRs) sim¬ 
ulations have served as a testbed for cosmological mod¬ 
els through comparisons with the observed large-scale 
distribution of galaxies. The first Horizon Run (HRl) 
was performed in 2008 and published in 2009 (Kim 
et al. 2009). The simulation box size was Lbox = 
6592 /i“*^Mpc, and the number of evolved particles was 
Np = 4120'^. The initial power spectrum was calculated 
by a fitting function from Einsenstein & Hu (1998), 
adopting a standard A cold dark matter (ACDM) cos¬ 
mology in a concordance with WMAP 5-year obser¬ 
vations (Dunkley et al. 2009). It produced eight non¬ 
overlapping all-sky lightcone data of halos and subhalos 
up to z = 0.6. We studied the non-linear gravitational 
effects on the baryonic acoustic oscillation (BAO) peak 
by measuring the changes in the peak position and am¬ 
plitude. In 2011, we performed even bigger simulations 
called Horizon Run 2 and 3 (HR2 and HR3, respec¬ 
tively; Kim et al. 2011). By adopting the same cos¬ 
mological model as in HRl, the initial power spectra of 
HR2 and HR3 were generated from the CAMB package 
(Lewis et al. 2000). The simulated galaxy distributions 
have been extensively exploited to measure both the ex¬ 
pected distribution of the largest structures for testing 
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cosmic homogeneity (Park et al. 2012) and the cosmic 
topology for constraining the non-linear gravitational 
effect on the halo density field (Choi et al. 2013; Kim 
et al. 2014; Speare et al. 2015). 

All the previous HRs have mean particle separations 
larger than 1 /i“*^Mpc, which has been sufficient for 
many cosmological tests. With much success in quanti¬ 
fying the non-linear gravitational effects on large-scale 
structures, recently we extended our research focus to 
galaxy formation studies. To model galaxies in simula¬ 
tions, we employed the subhalo-galaxy one-to-one cor¬ 
respondence model and abundance matching between 
subhalo mass function and the observed galaxy lu¬ 
minosity or stellar mass function (Kim et al. 2008). 
Most characteristics of observed galaxy distributions 
(in terms of luminosity functions and one-point den¬ 
sity distributions) are well reproduced by the model 
while observed abundances of galaxy clusters are not 
properly recovered from subhalos. The underpopula¬ 
tion of simulated galaxy clusters may come from the 
inefficient subhalo findings in cluster regions (Muldrew 
et al. 2011; for the various subhalo finding comparisons 
see Onions et al. 2012) or from the spatial decoupling 
between subhalos and galaxies. The latter may survive 
the tidal disruption longer due to more compact sizes 
through baryonic dissipation (Weinberg et al. 2008). 

In the ACDM cosmology, dark matter halos form hi¬ 
erarchically through the merger of smaller structures. 
These merger events can trigger star formation and 
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drive galaxy formation and evolution (Kauffmann et al. 
2004; Blanton & Berlind 2007). The merger history of 
galaxies has extensively been studied in semi-analytic 
models (SAMs; Cole et al. 1994; Kauffmann et al. 1997; 
De Lucia et al. 2004; Baugh 2006; Lee et al. 2014) for 
the last two decades. In SAMs, the gas heating and 
cooling rate are tabulated, and the resulting star forma¬ 
tion and supernovae feedback effects are implemented 
with some parametric prescriptions. Those parameters 
are fine-tuned to reproduce the correlation functions 
and/or luminosity functions of observed galaxies. Even 
though SAMs have achieved a great success in repro¬ 
ducing some observables, they require the introduction 
of a large number of parameters that are not necessarily 
physical. 

Another well-known empirical galaxy model, the halo 
occupation distribution (HOD; Berlind & Weinberg 
2002; Zheng et al. 2005, 2009) modeling, has been 
adopted to match the inner part of the galaxy correla¬ 
tion, which is attributed to satellite pairs inside a viri- 
alized halo. To distribute satellite galaxies in a halo, 
they empirically measure he probability number dis¬ 
tribution of satellites from observations. The HOD is 
simpler than SAMs, and widely used for the comparison 
between observed galaxies and simulated halos. 

The galaxy-subhalo correspondence (or the abun¬ 
dance matching; Kim et al. 2008; Trujillo-Gomez et al. 
2011; Rodriguez-Puebla et al. 2013; Reddick et al. 2013; 
Klypin et al. 2015) model is positioned between the two 
aforementioned models. It is much simpler than SAMs 
but based on more physical processes than the HOD. It 
originally models the satellite galaxy distribution from 
subhalo catalogs. Satellite galaxies in a galaxy clus¬ 
ter originally formed in situ isolated and merged into 
the cluster afterward. While falling into the potential 
well of the cluster, they experience a drag force by dy¬ 
namical friction (Zhao 2004), and they inevitably show 
spiraling inward orbital motions. After a certain time, 
they finally merge into the central galaxy. 

Although it seems reasonable to assume the presence 
of a satellite galaxy inside a subhalo as long as there 
is no galaxy-halo decoupling, it has been noted that 
some satellite galaxies may not have a host subhalo 
(Gao et al. 2004; Guo et al. 2011; Guo & White 2014; 
Wang et al. 2014). This could be tested by extensive 
hydrodynamical simulations to investigate the segrega¬ 
tion between satellite galaxies and their dark matter 
host (Weinberg et al. 2008). However, hydrodynamical 
simulations are expensive to run and still require much 
effort to reduce ambiguities in astrophysical processes 
and numerical artifacts. On the other hand, Hong et al. 
(2015) recently showed that if most bound particles are 
used instead of subhalos in the modeling, it is possible 
to identify such satellites without hosting subhalo. 

Therefore, we performed a new simulation in our se¬ 
ries, the Horizon Run 4 (HR4). This simulation, with 
improved spatial and mass resolutions with respect to 
the previous runs, retains a large number of particles. 
It is well-suited to study galaxy formation by producing 
merger trees. 


The outline of the paper is as follows. In Section 2 
and 3, we describe the simulation specifics and outputs 
of HR4, respectively. Mass function, shape and spin of 
virialized halos are dealt with in Section 4. The anal¬ 
ysis of two-point correlation functions and mass accre¬ 
tion history are given in Section 5 and 6, respectively. 
Summary and discussions are following in Section 7. 

2. GOTPM Code and Simulation 

2.1. Initial Conditions & Parallelism 

The simulation was run with an improved version of 
the GOTPM code (Dubinski et al. 2004). The input 
power spectrum is calculated by the GAMB package, 
and the initial positions and velocities of the particles 
are calculated by applying the second-order Lagrangian 
perturbation theory (2LPT) method proposed by Jenk¬ 
ins (2010). The gravitational force is evaluated through 
splitting the Newtonian force law into long- and short- 
range forces (for the Newtonian and Relativistic rela¬ 
tions, see Rigopoulos & Valkenburg 2015; Hwang et al. 
2012). The long-range forces are calculated from the 
Poisson equation in Fourier space for the density mesh 
built by the Particle-Mesh (PM) method. The short- 
range forces are measured with the Tree method. 

We parallelized the GOTPM code implementing MPI 
and OpenMP with a one-dimensional domain decom¬ 
position (z-directional slabs). We adopt a dynamic do¬ 
main decomposition, which sets the number of particle 
in each domain to be equal within one percent. Ac¬ 
cordingly, the slab width changes during the simula¬ 
tion run. By using a dynamic domain decomposition, 
one can easily identify the neighborhoods of a domain 
and establish communications between them. On the 
other hand, slab domains usually have greater surface- 
to-volume ratios than ordinary cubic domains (e.g., the 
orthogonal recursive bisection, Dubinski 1996), and so 
it has large communication size between domains. 

2.2. Non-recursive Oct-Sibling Tree 

We have employed a non-recursive oct-sibling tree 
(OST) for the tree-force update. The OST is a struc¬ 
tured tree of particles and nodes with mutual con¬ 
nections established by sibling and daughter pointers. 
Each particle has one sibling pointer, and each node 
has two pointers: one for its daughter and the other for 
the next sibling. 

First, we create a top-most node encompassing all 
particles. We define it as the zero tree level, and its sib¬ 
ling pointer is directed to the null value (Fig. 1). From 
the top-most node, we recursively divide each node into 
eight equal-sized cubic subnodes by increasing the tree 
level by one. If a sibling subnode contains more parti¬ 
cles than a predefined number, we divide the node fur¬ 
ther by increasing the tree level by one. The daughter 
pointer of the node is directed to the first sibling subn¬ 
ode, and the other subnodes are linked by their sibling 
pointers. If the node does not have any subnode, we 
make a chain of particles linked by their sibling point¬ 
ers, and we set the start and end of the chain connected 
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Figure 1. Example of the Oct-Sibling Tree structure. Boxes 
and circles represent nodes and particles, respectively. The 
black and blue arrows are daughter and sibling pointers, 
respectively. Each node has a daughter and sibling pointers 
while each particle has only a sibling pointer. 


to the previous and next sibling nodes (or possibly par¬ 
ticles), respectively. The last sibling at each local tree 
level is set to have its sibling pointer directed to the 
mother’s next sibling if it exits. If not, we also recur¬ 
sively climb the local tree until we find the next sibling 
of the current tree line. 

The advantage of the OST over the traditional oct- 
tree is a smaller number of pointers it employs and the 
needlessness of a recursive tree walk, which requires ad¬ 
ditional costs for such a stacking process to temporarily 
store information of the current recursive depth. The 
algorithm 2.1 is a pseudocode for the non-recursive tree 
walk with the OST. The tree walk is taken until the run¬ 
ning pointer, p, encounters the null value. During the 
tree walks the opening of a node is determined by the 
Open function. The tree-force update is done either by 
GroupForce or ParticleForce depending on the data 
type addressed by the pointer {p —>• type). These three 
functions play a pivotal role in tree walks. Open de¬ 
cides whether to go further into one deeper level (open¬ 
ing the node and going down to its daughter) or jump 
to the next sibling under the opening condition that 
9 > 9c, with 9c the predefined opening threshold. The 
GroupForce function calculates the gravitational force 
from the group of particles using the multipole expan¬ 
sion while ParticleForce calculates the gravitational 
force from the particle, p. Thanks to its cost efficiency, 
this kind of pseudocode is widely applied to our analy¬ 
sis tools such as a percolation method (Friend-of-Friend 
halo finding), peak findings in a Spherical Overdensity 
halo identification, and the two-point correlation mea¬ 
surement. 


2.3. Position Accuracy in GOTPM 

One of the key factors to determine the resolution of La- 
grangian codes is the spatial accuracy or, more specif¬ 
ically, the number of significant digits involved in the 


Algorithm 2.1: GotpmTreeWalk(p) 


while p ^ NULL 

"if (p —type) = NODE 

"if Open(p) = YES 


do < 


then {p := p ^ daughter 
, J call GroupForce(p) 
1 p := p —?► sibling 
J call ParticleForce (p) 

1 p := p — sibling 
comment: p is a running pointer. 


then 


else 


particle position. Usually, a single-precision floating¬ 
point type has been applied to save the position of a 
particle because the four-byte single precision is suffi¬ 
cient for small simulations. However, as the number 
of particles in simulations increases, the position accu¬ 
racy from single-precision begins to deteriorate. Since 
the roundoff error of a single precision variable A is 
eroundoff(-4) ~ 10“^^!, the maximum roundoff error of 
the single-precision position with respect to the mean 
particle separation is 


^roundoff 


10 


-7 




For example, if the total number of particles is 6300^ as 
in the HR4, the maximum roundoff position error lies 
at the level of a few sub-percent of the mean particle 
separation, or eroundoff(7' max/^mean) 10-3. 

On the other hand, in the HR4 as well as in the HR2 
& 3, we separate the position of a particle (r) into two 
vectors as 


r = L + d, 


( 2 ) 


where L and d are the Lagrangian position and dis¬ 
placement from the Lagrangian position of a particle, 
respectively. We set the particle index by a row-major 
order in the Lagrangian configuration and, therefore, it 
does not require additional memory space to compute 
L. Since the displacement of the simulated particle over 
the entire HR4 simulation run is less than ten times 
of the mean particle separation (dmax ^ lOdmean), the 
maximum position error in the HR4 is 


^roundoff 


10 


7 ^max 


10 


-6 


( 3 ) 


In this way, we significantly enhanced the accuracy of 
the particle position without using any additional mem¬ 
ory space. 


2.4. Simulation Specifics 

The HR4 was performed on the supercomputer of 
Tachyon-II installed at KISTI (Korea Institute of Sci¬ 
ence and Technology Information). We used 8,000 CPU 
cores over 50 straight days from late November in 2013 
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to early February in 2014. Even with several system 
glitches over the allocated time period, we succeeded 
to complete the simulation in about 50 days for the 
gravitational evolution of 6300^ particles in a periodic 
cubic box of a side length Lbox = 3150 /i“^Mpc. The 
starting redshift is Zi = 100, which is chosen for par¬ 
ticles not to overshoot one grid cell spacing (Lukic et 
al. 2007) in setting the initial conditions. This high 
initial redshift, combined with 2LPT, ensures an accu¬ 
rate power spectrum and mass function measurement 
at z = 0 (L’Huillier et al. 2014). The simulation took 
2000 steps to reach the final epoch of Zq = 0. The mean 
particle separation is set to dmean = 0.5 /i“^Mpc and 
the corresponding force resolution is O.ldmean- 

We adopted a standard ACDM cosmology in concor¬ 
dance with WMAP 5-year. This choice of cosmology 
was made for consistency with various observations in¬ 
cluding SDSS as well as the previous HRs. Specifically, 
the matter, baryonic matter, and dark energy densities 
are = 0.26, ^lb,o = 0.044, and = 0.74, respec¬ 
tively. The current Hubble expansion is Hq = 100 h 
km/s/Mpc, where h = 0.72. The amplitude of the ini¬ 
tial matter perturbations is scaled for an input bias 
factor, bs = l/cs = 1-26, where 

^1 = ^ j k^P{k)\W{kR^)\'^dk, (4) 

and i?8 = 8/i“^Mpc. Here, we used the spherical top- 
hat filter W{x) = 3(a:sina: — cosx)/a:^ in fc-space. The 
particle mass is nip ~ 9 x 10® and the mini¬ 

mum mass of halos with 30 member particles is about 
Ms ~ 2.7 X lO^i 

Figure 2 shows the evolution of the non-linear mat¬ 
ter power spectrum obtained during the simulation run 
at several redshifts. The dotted lines are the expected 
linear power spectra, while the solid lines are the sim¬ 
ulated matter power spectra at the same redshift. The 
typical non-linear evolution effect can easily be seen on 
small scales, where the amplitude of the power spec¬ 
trum is greater than the linear prediction due to the 
gravitational clustering. 

Figure 3 shows a part of the density map of the HR4 
at z = 0, where a cluster develops at the center through 
the mergers of several neighboring overdensity clumps. 
One may clearly see void regions (painted in dark blue) 
with a size of a few tens of /i“^Mpc. Some overdense 
blobs are embedded in the connection of multiple fila¬ 
mentary structures. 


3. Outputs 

In this section, we describe the main products of the 
HR4 simulation. They are available at http://sdss. 
kias.re.kr/astro/Horizon-Run4/. 

3.1. Snapshot and Past Lightcone Space Data 

We have saved snapshot data of particles at twelve red- 
shifts: z = 0, 0.05, 0.1, 0.15, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 
1, and 4. Each data set contains the particle position, 
velocity, and eight-byte integer ID index. 



k (h/Mpc) 


Figure 2. Matter power spectra from the HR4 simulation 
{solid) and from linear theory {dotted). Color: z = 100 
{black), 3.6 {blue), 0.32 {red), and 0 {magenta). 


To generate the past lightcone space data, we put 
an artificial observer at the origin {x, y, z) = (0,0,0) of 
the simulation box. At each time step, we calculate the 
comoving distance from the observer using 


dc 


-I 

Ho Jo 


E{z') 


dz'. 


(5) 


where 


E{z) = 

+ ^)^ + (1 ~ ^m,0 ~ ^A,o)(l + ^)® + IIa.O- 

( 6 ) 

Then, we search for particles located in a comoving 
shell, whose inner and outer boundaries at the i-ih step 
are dc,i — Adc,i/2 and dc,i + Adc,i+i/2, respectively, 
where Adc,i+i = {dc,i+i — dc,i)l2. We utilize the pe¬ 
riodic boundary conditions by copying the simulation 
box to extend the all-sky past lightcone space data up 
to r = 3150 /i“^Mpc, which corresponds to z ~ 1.5. 

Due to the finite step size, several undesirable events 
may be encountered. If a particle crosses the shell 
boundary between two neighboring time steps, it can 
be missed or be counted twice in the lightcone space 
data. Therefore, we set a buffer zone laid upon both 
sides of the shell. The width of the buffer zone is deter¬ 
mined to be equal to the maximum displacement taken 
by a particle in a time step. Using these buffer zones, 
we can catch these crossing events. A crossing particle 
can simultaneously be detected in two contacting shells 
or two adjacent buffer zones, and we simply merge the 
duplicated particle by averaging position and velocity 
in the lightcone space data. 
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Figure 3. Simulation density slice map at z = 0. High-density regions are painted with bright color. The width of the slice is 
7 h“^Mpc. The two subfigures are arranged for cascaded zoom-in views of a cluster at the center of the box in the bottom 
part of the figure. We put a scaling bar on the bottom of each panel. 
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Figure 4. Bottom- Distribution of galaxies with A4r < —22.35 and 0.45 < 2 < 0.6 in the BOSS-CMASS volume-limited 
sample (Choi et al. 2015). Galaxies are selected from a strip with —2° < dec. < 2° and 130° < R.A. < 230° out of the entire 
BOSS survey area. Top-. The ‘galaxies’ in the mock BOSS-CMASS survey performed in the HR4 simulation. The absolute 
magnitude of the observed galaxies is calculated with the reference redshift of 2 = 0.55 and is used to produce the sample 
with a constant number density. The HR4 PSB subhalos are selected as galaxies in accordance with the galaxy-subhalo 
correspondence model, and the galaxy number density at a given redshift is matched with the observed one by adjusting 
the low-mass cutoff. 


In both the snapshot and past lightcone data, we ap¬ 
ply the Ordinary Parallel Friend-of-Friend (OPFOF), a 
parallel version of FoF code to identify virialized halos. 
The standard percolation length is simply adopted as 
^link = 0.2(ii„ean- The halo position and peculiar ve¬ 
locity are given as the average position and velocity 
of the member particles. Then we apply the physically 
self-bound (PSB) subhalo finding method (Kim & Park 
2006) to identify subhalos embedded in the FoF halo. 
It employs a negative total energy criterion and spher¬ 
ical tidal boundaries to discard particles from subhalo 
candidates. 

As a representative example, Figure 4 compares the 
volume-limited galaxy sample from BOSS-CMASS with 
r-band magnitude limit Air < —22.35 at 0.45 < 2 < 
0.6 {bottom) with the mock galaxy sample from the HR4 
built by the PSB subhalo-galaxy correspondence model 
{top-, Kim et al. 2008). 

3.2. Halo Merger Data 

To build the merger trees we detect halos and subhalos 
at 75 equally-spaced sparse time steps from 2 = 12 to 
0. The step size is set to be comparable to the rota¬ 


tional period (i.e. dynamical timescale) of Milky-Way- 
size galaxies. Halo merger trees are then built by trac¬ 
ing the gravitationally most bound member particles 
(MBPs) of halos. If a halo does not contain any for¬ 
mer MBP, we select a new MBP among the member 
particles of the halo. If one MBP is found, we assume 
the halo to be a direct descendant of the halo. If the 
halo hosts multiple former MBPs, we treat the halo as 
a merger remnant and those ancestor MBPs (or halos) 
are linked to the remnant creating a halo merger tree. 
These merger trees will be extensively used to build 
mock galaxies and to compare with observations (Hong 
et al. 2015). Of course, due to the halo mass resolution 
of the HR4 {Mg = 2.7 x 10^^ /i^^Mq), we are unable 
to resolve mergers of sub-Milky-Way-mass (sub)halos. 

4. Properties oe FoF Halos 
4.1. Multiplicity Function 

The multiplicity function is defined as 

f( j — ^ dn(M, 2 ) 

P 6 ( 2 ) dln[i/(T(M, 2 )]’ 


( 7 ) 
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ln(l/(7) 


Figure 5. Bottom-. Multiplicity functions of HR4 FoF halos 
at « = 0, 0.5, 1, 1.99, 4, and 5. The solid curve is pro¬ 
posed by Bhattacharya et al. (2011), Bll. Top: Fractional 
deviations of the simulated {symbols) and modeled {lines) 
multiplicity functions from Bll. 


where n{M, z) is the cumulative halo mass function at 
z, Pb{z) is the background matter density, and a{M,z) 
is the density fluctuation measured on the mass scale 
of M. For a given power spectrum P{k), the density 
fluctuation is estimated as 

a\M,z) = I k^P{k)\W{kR{M,z))\^dk, (8) 


where 


R{M, z) 


/ 3M \ 

\iTTpb{z)) 


(9) 



z 


Figure 6. Redshift-dependent y- and amplitude-corrections 
in Equation 11 showing the best fit to the HR4 FoF mul¬ 
tiplicity functions. The dotted lines are the analytic fitting 
functions as shown in Equations (12) and (13). 


multiplicity functions is substantial (Lukic et al. 2007), 
and the overall amplitude seems to increase as the red- 
shift decreases. 

We fit the simulated multiplicity function with a vari¬ 
ant of the Bll function with an amplitude changing 
with redshift as 

/Kim(XL,2) = ip{z)fBll{Xh{M,z) - Xb{z)). (H) 


and Di{z) is the growing mode of the linear growth 
factors computed as 

D,{Z) = l^rn,0E{z) (10) 

Figure 5 shows the multiplicity function from the 
HR4 as well as a number of previous fitting models 
of the multiplicity function (see Table 1 and references 
therein). Top panel shows the deviations of multiplicity 
functions with respect to the model described in Bhat¬ 
tacharya et al. (2011), hereafter Bll. For clarity, in the 
cases of Crocce et al. (2010) and Manera et al. (2010) 
we only show the fitting function obtained at z = 0. All 
the previous fitting functions significantly deviate from 
each other at high mass scales. This may be produced 
by the exponential cut off producing large noises in fit¬ 
ting the steep slope. From the simulated multiplicity 
function, we can clearly see the redshift change, and 
therefore a single functional form may not be sufficient. 
For large values of ln(l/cr), the redshift evolution of 


Here xl{M, z) = y/q6c/cr{M, z), where (5c = 1.686 is the 
density contrast at the collapse epoch in an Einstein- 
de Sitter universe, g is a fitting parameter in the Bll 
function (see Table 1), and Xs{z) and (p(z) are redshift- 
dependent y- and amplitude-corrections, respectively. 
The value of 6c in an Einstein-de Sitter Universe is 
1.686, and slightly depends on the cosmology. How¬ 
ever, for consistency with previous work, we will use 
this constant value of 1.686 (e.g., Bhattacharya et al. 
2011). We fit the simulated multiplicity function with 
the least-square minimization and obtain the empirical 
fitting function as 

Xs(z) = 0.09 tanh^(0.9z)-b 0.01 (12) 

ip{z) = exp (-^) + 0.025. (13) 

Figure 6 shows the redshift evolution of Xs{z) (left) 
and (/3(z), respectively. At high redshifts, the HR4 sim¬ 
ulation underpopulates halos compared to Bll, while 
the HR4 has more halos than Bll in the recent epoch. 
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Table 1 

Description of fitting models of the multiplicity function 


Model 

/(T2) 



Parameters 

Sheth & Tormen (1999)* 

(1 + X exp 

1 - 

1 

{A,p,q) = (0.3222,0.3,0.707,0.3) 

Jenkins et al. (2001) 

A exp liner 


(A, a, 6) = (0.315,0.61,3.8) 

Warren et al. (2006) 

A (cr~“ + 6) exp 

C ■ 


(A, a, b, c) = (0.7234,1.625, 0.2538,1.1982) 

Tinker et al. (2008)^ 

A (cr~“ + 6) exp 

. 

C 


(A, a, b, c) = (0.745,1.47,0.250,1.19) 

Crocce et al. (2010)^ 

A (cr~“ -I- 6) exp 

c 


(A, a, b, c) = (0.58,1.37, 0.30,1.036) 

Manera et al. (2010)*^ 

(1 + X exp 

1 - 1 

1 

{A,p,q) = (0.3222,0.248,0.709) 

Bhattacharya et al. (2011)* 

A'y^x’"(l + X“^^)exp 

r x^ 

2 

{A,p,q,r) = (0.333,0.807,0.788,1.795) 

Angulo et al. (2012) 

/ B \ 9 

A f — + 1 j exp 

C 



(A, B, C, q) = (0.201, 2.08,1.172,1.7) 

Watson et al. (2013) 

A (cr~“ -I- 6) exp 

-a 


{A,a,b,c) = (0.589,2.163,0.479,1.210) 


* X = where Sc = 1.686 is the density contrast at the collapse epoch in an Einstein-de Sitter universe. 

1 Only the case at 2 = 0 is given here. 



ln(l/CT) 


Figure 7. Simulated and modeled multiplicity functions with 
respect to our new fitting model (/Kim(cr, 2)). Same symbol 
and color conventions as in Figure 5. 

Also, as we move to higher redshift, Xs{z) increases and 
reaches about 0.098 at very high redshift. 

Figure 7 shows the scatter of the simulated multi¬ 
plicity function (symbols) over our htting model (/Kim) 
with various former models (lines) for comparison. 
While other fitting models match the simulated multi¬ 
plicity function only at small scales (ln(l/cr) < 0.3), our 
new fitting model agrees with the HR4 on most scales 
in the redshift range between z = 5 and 0, within a 
2.5% fluctuation level. 

The origin of the redshift evolution of the multiplic- 



Figure 8. Shape distributions of FoF halos in various mass 
samples shown in the q-s diagram at 2 = 0. We mark iso¬ 
density contours enclosing 25% halos around peak distribu¬ 
tions. 
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ity function is not clearly known but it might be partly 
explained by the following arguments. First, even the 
2LPT may be somewhat insufficient to generate accu¬ 
rate initial conditions of the simulation. Such effect 
would be avoided by introducing higher-order approx¬ 
imations (for comparison between the Zel’dovich ap¬ 
proximation and 2LPT, see Crocce et al. 2006). How¬ 
ever, the target redshifts to measure the halo mass func¬ 
tion are sufficiently lower (z < 5; see the discussions 
made by Tatekawa & Mizuno 2007) compared to the 
initial epoch of Zi = 100 and, consequently, the red- 
shift evolution may not be caused by numerical tran¬ 
sients. Second, it may be due to the limitation of the 
linear perturbation theory or the linear growth model 
in calculating a{M,z) after the nonlinear gravitational 
clustering begins to enter. The multiplicity function as¬ 
sumes that there is no other redshift dependence than 
the matter fluctuation, a{M,z) = D{z)a{M). But as 
the non-linear growth becomes significant at lower red- 
shifts, one should consider the effect of non-linear grav¬ 
itational evolution of density fields. 


4.2. Halo Shape 

4.2.1. Structure 


The shape tensor Sij of an FoF halo is defined 

Nm 

Sij = -Xj), (14) 

k 


where i and j are structural axes, x is the position of 
the center of the halo mass and Nm is the number of 
member particles. The three eigenvalues of the shape 
tensor < r 2 < ri are respectively the lengths of the 
minor, intermediate, and major axes of the correspond¬ 
ing ellipsoid. The prolateness and sphericity of a halo 
are defined as 


r2 



s = — 
ri' 


(15) 

(16) 


A halo is respectively defined as prolate, oblate, or 
spherical if 


g < 1, 

(17) 

g ~ 1, s 1, 

(18) 

s ~ 1. 

(19) 


Figure 8 shows the probability distributions of {q, s) 
with a contour containing 25% of halos around the peak 
distribution in four different mass samples of FoF halos. 
For halo samples more massive than 2 x 10^^ h~^MQ, 
we can clearly see that halos become more prolate as 
the mass increases, in agreement with theoretical pre¬ 
dictions (e.g. Rossi et al. 2011). Less massive halos 
with Mg <M<5xlO^^ h~^MQ have their distribu¬ 
tion substantially shifted to the lower-right corner in 
this diagram, i.e., are more oblate than more massive 
samples. This may be fully explained by the particle 
discreteness effect, which will be described in the next 
section. 



Figure 9. Resolution dependence of the roundness parameter 
at z = 0 for HR4 (blue), GRl {green), and GR2 {magenta) 
containing 25% of halos around peak probabilities. Bot¬ 
tom: Roundness parameter as a function of halo mass. A 
vertical bar marks the mass scale equivalent to the mass of 
1000 particles (lO^rrip) combined. The best-fitting functions 
of TZ{M) from halos with lower-mass cutoff lO^rrip (7?.fooo; 
yellow) and 5 x lO^rrip (77.5ooo; red dash) are shown. Top: 
Deviation of the roundness parameter from 7?.5ooo with re¬ 
spect to the halo mass. Here we do not show the distribution 
below the mass of lO^rrip. 
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4.2.2. Roundness 

We now investigate the FoF halo shape from a different 
angle. First, we define the roundness as 


( 20 ) 

To measure the resolution effects on halo shape, we 
ran two additional higher-resolution simulations called 
Galaxy Run 1 (GRl) and Galaxy Run 2 (GR2). These 
simulations used 2048^ particles. We employed the 
same cosmological model but different simulation box 
sizes = 512 h“^Mpc & = 256 h~^Mpc). 

The mean particle separations of GRl and GR2 are 
c^mean = 0.25 and 0.125 /i“^Mpc, respectively, while the 
corresponding force resolutions are changed to 0.025 
and 0.0125 /i“^Mpc accordingly. Therefore, the mass 
and force resolutions are quite enhanced with respect 
to the HR4. 

Figure 9 shows the distribution of TZ of FoF halos at 
z = 0. In each simulation, at scales of M ^ lO^Wp, TZ 
tends to be independent of the simulation resolution. 
On the other hand, at small mass scales (M < lO^rrip) 
each simulation seems to underestimate TZ, probably 
due to the small number of particles. Hoffmann et al. 
(2014) examined the discreteness effect on a modeled 
halo for a given shape and found that the required 
number of particles should not be less than 1000 for 
a reliable shape determination (see Fig. A2 in their 
paper). 

From our three simulations we found a fitting formula 
of the roundness parameter as a function of the halo 
mass for massive halos with M > lO^Wp: 


^1000 logj^Q 


M 


107 h-iM0 


X exp 


M 


^Klogio 


0/ J 


( 21 ) 


where = (0.55,0.28) (yellow line in Figure 9). 

One should note that this fitting is only valid for M > 
1.4 X 10^7 which is set by the combined mass 

of 1000 particles of the GR2. We do not find any turn¬ 
around mass scale of TZ in the available mass range in 
this study. If we only consider halos with M > 5 x 
lO^rUp, the distribution of TZ follows 


”^5000 (-^) — Tf^logio ^2^012 /J-IM 0 ) (22) 

where {A-n, B-n) = (—0.07,0.68) (red dashed line in 
Figure 9). Similar to the case of TT-foooi the above fit¬ 
ting is only valid for M > 7 x lO^^ /i“iM 0 . Both TZ^qqq 
and 77 .|qoq describe well the change of TZ, but they di¬ 
verge below M = 7 X 10^^ /i“ 7 ]V[ 0 , the mass scale of 
about 5 X lO^TOp of GR2. Therefore, we may need a 
simulation with a higher mass resolution than GR2 to 
see which fitting formula describes the roundness pa¬ 
rameter of low-mass halos. 



Figure 10. Change of peak position in the TZ distribution at 
several redshifts in the HR4. The lower bound of the x-axis 
corresponds to 10 ^Trip. 


Figure 10 shows the redshift evolution of TZ in the 
HR4. At higher redshift, halos tend to have a smaller 
value of TZ for a given mass. However, it is important to 
note that this tendency does not guarantee a possible 
shape evolution of virialized halos because halos also 
grow in mass with time. 


4.3. Halo Orientations 

In this section, we study the angle between the halo 
rotational and structural axes. The directional angle 
between them is calculated as 

di = cos“^ |J • n|, (23) 


where J is the normalized rotational axis and is the 
unit vector in the direction of the structural axis i. We 
define the probability distribution function of the direc¬ 
tional angles 


p'{e) ^ 


dP( 6 ») 
d cos 9 ’ 


(24) 


where P{9) is the cumulative probability of a direc¬ 
tional angle greater than 9. Then, for a random orien¬ 
tation p'(9) is uniform over the angle of 0° < d < 90°. 

Figure 11 shows the relations between the rotation 
and halo axes as a function of halo mass. The rotational 
axis tends to be orthogonal to the major axis (bottom 
panel), which means that halos tend to swing around 
their major axis. Moreover, from the upper two panels, 
it can be noted that the rotational axis is more aligned 
with the minor axis than the intermediate axis. This 
alignment becomes stronger as the halo mass increases. 
In addition, we find that this tendency still holds for 
low-mass halos with M <10^mp, implying that the halo 
rotation is less affected by the mass resolution limit 
than the halo shape. 













10 


J. Kim et al. 


0) 

0) 

s- 

M 

0) 

T3 


80 - 

60 - 75 % 


40 
20 - 




50X 


0 Tf- fl tlr+p^ 

80 - 


Bx • pa-ij*®#!*! 


.i'l'il !l 


JIVf 


[nil 



75% 




WMfcs' , 


75% 


10>2 


10>3 IQX 

M(h-iMg) 


10>6 


Figure 11. Orientations of the rotational axis with respect to 
the major {bottom), intermediate (middle), and minor (top 
panel) axes at 2 = 0. Probability contours are drawn around 
the peak position at each mass bin enclosing 25%, 50%, 
and 75% of halos respectively. Most halos are positioned at 
01 ~ 90°, 02 ~ 0°, and 03 ~ 0°. 


5. Two-Point Correlation Function 


In this section, we implement the effects of redshift- 
space distortions on the clustering of mock galaxies and 
measure the change of clustering in the radial (tt) and 
tangential (cr) directions. In the 3-dimensional space, 
the radial separation between two points (ri and r 2 ) is 
defined as 


TT = 


\di2 ■ i?i2| 

\Ri2\ 


(25) 


where R 12 = (ri -|- r2)/2 and <ii 2 = Ti — r 2 . The 
tangential distance between them is simply obtained 
with 


cr = a/ di2 ■ di2 - TT^. (26) 


The correlation function of a point set can easily be 
calculated by Hamilton’s method (Hamilton 1993): 


DD{a,Tr)RR{a,Tr) 
DR{a, 7r)2 


(27) 


Here, DD is the number of pairs of real points, DR is 
the number of cross pairs between the real and random 
points, and RR is the number of pairs of random points 
at the two-dimensional separations of a and tt. 

We use the PSB subhalo catalog from the HR4 snap¬ 
shot at z = 0 as our mock galaxy sample. By adopt¬ 
ing the far-field approximation and using the periodic 
boundary condition of the HR4 simulation, we produce 
the redshift-space distortion in the ^-direction. 


x' = x+^, (28) 

rio 

where Hq is the Hubble parameter at z = 0 and Vx is 
the peculiar velocity along the x-axis. Since we adopt 
the far-field approximation, tt is the position difference 
in the a;-axis and a is the separation in the y-z plane. 
We then construct a mass-limited mock galaxy sample 
with PSB subhalos satisfying M > 2.60 x 10^^ 

The average number density of the mass-limited PSB 
subhalo sample is n = 1.48 x 10“^ /i^Mpc“^, which 
is comparable to the number density of the volume- 
limited sample of the SDSS Main galaxies with absolute 
magnitude limit of Mr — Slogj^p^ < ~21 (Choi et al. 
2010 ). 

Figure 12 shows the effects of redshift-space distor¬ 
tions on the correlation map. The left panel shows 
the correlation of our mock galaxy sample in real space 
while the effects of redshift-space distortion are applied 
in the right panel. The shape of ^{TT,a) is distorted 
along the line-of-sight (LoS, or in the 7r-direction). At 
the very center, the finger-of-god effect can be seen as 
spikes stretching along the 7r-direction (for a better view 
around the center, see Figure 13). On the other hand, 
on larger scales, the correlation function along the LoS 
contracts to the smaller scale. 

The position of the BAO peak in real space can be 
estimated from the linear correlation function 


r{r) = 


27r^ 


'k^P(kf^dk^ (29) 


kr 







































Horizon Run 4 


11 







0 


-1 




-2 


-100 -50 0 50 

o (h ’ Mpc) 




-100 -50 0 50 100 

a (h ' Mpc) 


(a) real space 


(b) redshift space 


Figure 12. Correlation functions of mock galaxies measured without {left) and with {right) redshift-space distortion effects. 
The radius of each circular region is 130 /i“^Mpc, and the solid circle marks the BAO peak position (upeak — 107 /i“^Mpc). 
The color bar marks the correlation in logarithmic spacing. 
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Figure 13. Same as the right panel of Figure 12, but zoomed 
in to clearly show the finger-of-god effect. 


For the WMAP 5-year standard ACDM cosmology, 
the BAO peak in real space is located at Upeak — 
107 h“^Mpc, shown as a solid circle in Figure 12. 

Figure 14 shows the two-point correlation functions 
for different values of the directional cosine to the LoS 
direction /r in real space {top) and redshift space {bot¬ 
tom panel). In real space, the correlation function 
around the BAO peak is independent of the directional 
angle {6) because of the isotropic distribution. On the 
other hand, the correlation functions measured in red¬ 
shift space are increased as 0 increases, because galaxy 
pairs are stretched along the LoS. It is worth to note 
that the BAO peak in the tangential direction {6 = 90°) 
cannot be detected. Moreover, the correlation functions 
along the LoS has a peak with a height nearly zero while 
correlation functions for 6 < 30° are less than zero on 
scales below the peak position. 

Figure 15 shows the average correlation function over 
the directional cosine, 

(30) 

Jo 

In both real and redshift spaces, the BAO peak from 
the HR4 is broadened and shifted toward small scales 
compared to a simple estimation of the linear correla¬ 
tion function of biased objects 

Clinear,bias(^) — b ^linear(^); (^1) 

where b = 1.14 is the bias factor. This is due to the 
nonlinear gravitational evolution of galaxies. In redshift 
space, the BAO peak is further broadened, and it is 
hard to clearly find the position of the BAO peak. 

6. Mass Accretion History 

We use merger trees to study the mass accretion history 
of halos in several mass samples. We dehne the mass 
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60 80 10 ® 1 . 2 x 102 

r (h“* Mpc) 


Figure 14. Correlation functions between PSB halo pairs 
separated along the directions of 0 = 10, 20, 30, 40, 50, 
60, 70, 80, and 90 degrees in real space {top) and redshift- 
distorted space {bottom panel). In the bottom panel, the 
top-most line is the correlation oi 9 = 80°, and the corre¬ 
lation function increases as 6 decreases. The dotted line is 
the linear prediction with a bias factor 6 = 1.14. As the 
direction cosine (/r = cos 6) increases, the noise of the cor¬ 
relation functions is decreasing because the number of pairs 
along the given direction increases (A^pair oc /r ). The cor¬ 
relation functions along 9 = 90° are not shown due to big 
noise. 



60 80 102 1 . 2 x 102 

r (h“' Mpc) 


Figure 15. Correlation function averaged over the directional 
cosine. The thick solid line is the averaged value of the 
correlation functions in real space {top) and redshift space 
{bottom panel). The dotted line is the linear prediction with 
bias b = 1.14. 

accumulation history as 

(32) 

Mq 

where Mq is the final halo mass at z = 0. The 
half-mass epoch {zi/ 2 ) is defined as the time when 
^{zi/ 2 ) = We measure the evolution of the 

halo mass along the major descendant trees and show 
the results in Figure 16. It can be seen that the half¬ 
mass redshift tends to decrease as the final halo mass 
increases. For example, low-mass halos with 10^^ < 
Mo//i“^M 0 < 3 X 10^^ tend to have their half-mass 
around Zi /2 — 1 on average, while more massive halos 
with 10^"^ < Mo/h“^M 0 < 5 x 10^“^ tend to have a 
later half-mass epoch Zi /2 — 0.5. This result is consis¬ 
tent with the observations that galaxy clusters formed 
relatively recently (in terms of the epoch when a cluster 
obtains half of the current mass) while individual satel¬ 
lite galaxies seem to form at relatively higher redshift. 

We empirically fit the log-linear function of redshift 
to T'(z) as 


T'(Mo, z) = exp [-i/i(Mo)z], (33) 

and found a best-fit relation as 

V'(Mo) 0.32 log^o 

This fitting function reproduces the distribution for 
10^^ < M//i“^M 0 < 5 X 10^® quite well at an early 














Horizon Run 4 


13 



0 18 3 4 5 

Z 


Figure 16. Evolution of halo mass with redshift for several 
mass samples. In the top panel, we show the change of ip 
with redshift, while is shown in the bottom panel. Lines 
and shaded regions mark the mean and la distributions of 
mass history. For each sample, we cut the data below the 
mass resolution of the simulation. 


epoch (z > 0.7). On the other hand, at low redshifts 
(z < 0.7) the accelerated expansion driven by dark en¬ 
ergy begins to overpower the gravitational attraction 
and, therefore, halo mergers are suppressed. The sharp 
increase of ip near the current epoch is caused by a 
numerical noise. Note that Dekel et al. (2013) also 
found an exponential form of mass growth, although 
their mass accretion rate ('0(10^^ //.“^M©) = 0.76) is 
slightly higher than ours (0(10^^ = 0.56). 

We want to point out that the specific mass accretion 
rate per unit redshift interval, defined as 


dM 


Mo d^i 

Mdz 


M dz 


0 (Mo), 


(35) 


stellar mass evolution of a halo is fully determined by 
the evolution of its total mass. In this case, the spectral 
indices of stellar mass-to-total mass, stellar mass accre¬ 
tion rate-to-stellar mass, and star formation efficiency- 
to-total mass, respectively defined as 


l{Mo,z) 

dlnM* 

dlnM 

(39) 

l3iMo,z) 

_ dlnT* 

“ dlnM* 

(40) 

e(Mo,z) 

din 5* 

“ dlnM’ 

(41) 


are fully determined by the redshift and the final halo 
mass. By applying the galaxy-subhalo correspondence 
model to relate between halo mass and galaxy lumi¬ 
nosity, Kim et al. (2008) showed that stellar luminos¬ 
ity (or stellar mass if a constant M*/L* is assumed) 
shows a good relation to the halo mass with a power- 
law index 7 ^ 0.5 for the SDSS main galaxy sample 
when M > 5 X lO^^/i^^M©. A similar slope was re¬ 
ported by Kravtsov et al. (2014) from the BCG sam¬ 
ples. On the other hand, Abramson et al. (2014) re¬ 
ported /3 ~ —0.3 for SDSS DR7 galaxies with 9.5 < 
logio(M*// i-^Mq) < 11.5. 

The spectral index of the stellar mass accretion rate- 
to-total mass can be expressed as a combination of the 
above spectral indices: 


r]{Mo,z) 


dlnT„ 

dlnM 

l3{Mo,zh{Mo,z) 


;(Mo, z) + 


ip{Mo) 


dlnA(z) 

Tz 



(42) 

(43) 
,(44) 


where E{z) was defined in Equation ( 6 ). The effect of 
the parameters on the relative star formation efficiency 
is shown in Figure 17. As can be seen from the figure, 
the relative star formation efficiency is higher, or the 77 
is getting smaller for more massive halos. 


7. Summary 


is roughly constant with redshift and depends only on 
the current sample mass. Then the specific mass ac¬ 
cretion rate per unit physical time can be calculated 
as 

Tm{Mq, z) = 


We now introduce a star formation efficiency, which 
is defined as the ratio of mass accretion rates between 
the stellar and total masses of halos as 

K{Mo,z) = ^, (38) 

I M 

where T* = dM-^/M^,dt is the specific stellar mass ac¬ 
cretion rate. As a simple case, we assume that the 


dM 
Mdt 
0(Mo) E{z) 
Hq 1 z 


(36) 

(37) 


We ran a new cosmological N-body simulation called 
the Horizon Run 4 (HR4) simulation. By adopt¬ 
ing a standard ACDM cosmology in concordance with 
WMAP 5-year observations, the HR4 simulates a pe¬ 
riodic cubic box of a side length, Lbox = 3150/i“^Mpc 
with 6300^ particles. With its wide range of mass and 
length scales, the HR4 can provide the cosmology com¬ 
munity with a competitive data set for the study of 
cosmological models and galaxy formation in the con¬ 
text of large-scale environments. 

The main products of the HR4 are as follows. First, 
we saved the snapshot data of the particles within 
the whole simulation box at 12 different redshifts from 
z = 4 to 0. We also built a past lightcone space data 
of particles that covers the all-sky up to z ~ 1.5. They 
can be used to study the evolution of the gravitational 
potential and the genus topology as well as large-scale 
weak lensing analysis. Moreover, we constructed the 
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Figure 17. Relative star formation efficiency scaled with the 
current efficiency, &*o. Clockwise from the lower-left panel, 
the halo masses are Mq = Mo = 10 ^^ 

Mo = 10 ^® h“^M 0 , and Mo = 10 ^® h“'^M 0 , respectively. 
In the legend, we list the values of rj from the bottom curve. 


merger trees of Friend-of-Friend halos from 2 ; = 12 to 
0 with their gravitationally most bound member parti¬ 
cles. They can be used to study galaxy formation and 
bridge the gap between theoretical models and observed 
galaxy distributions. 

We tested the HR4 in various aspects, including the 
mass/shape/spin distributions of FoF halos, two-point 
correlation functions of physically self-bound subhalos, 
and mass evolution of FoF halos. The results of our 
test are summarized as follows: 

1. We found that the abundance of massive FoF ha¬ 
los in the HR4 is substantially different from var¬ 
ious fitting functions given in the previous litera¬ 
ture. We also found strong evidence for a redshift 
dependence of the mass function. We proposed 
a new fitting formula of the multiplicity function 
that reproduces the redshift changes of amplitude 
and shape of multiplicity functions within about 5 
% errors. 

2. We conhrmed the finding of previous studies that 
FoF halos tend to rotate around the minor axis. 

3. The two-point correlation function measured in 
real space is isotropic. However, due to the non¬ 
linear evolution of galaxies, the location of the 
baryonic acoustic oscillation peak is shifted toward 
smaller scale than the prediction from the linear 
correlation function. On the other hand, in red¬ 
shift space the BAO peak can be seen only in the 
two-point correlation function along the perpendic¬ 
ular direction, with a much-broadened width and 
increased height. We emphasized that it is im¬ 


portant to use massive simulation data to study 
the non-linear evolution of BAO features and the 
connection between observations and cosmological 
models. 

4. We found that more massive halos tend to have 
steeper mass histories, and the mass accretion rate 
per unit redshift is roughly constant during early 
epoch before dark energy domination. By adopting 
simple power-law models for the stellar mass and 
star formation efficiency, we found that massive ha¬ 
los tend to have a higher star formation efficiency. 

All aforementioned main products of the HR4 
are available at h.ttp://sdss.kias.re.kr/astro/ 
Horizon-Run4/. 
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